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Hole simulations, either with or without shift). The code is built upon an adjusted first- 
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In a recent paper Kiuchi and Shinkai have analyzed numerically the behavior of many 
'adjusted' versions of the BSSN system. This is a follow-up of a former proposal [2j for using the 
energy-momentum constraints to modify Numerical Relativity evolution formalisms. An important 
point was to put the constraint propagation system (subsidiary system) in a strongly hyperbolic 
form, so that constraint violations can propagate out of the computational domain. As a further 



step, there is also the possibility of introducing damping terms, which would attract the numerical 

o 

solution towards the constrained subspace. 

At the first sight, one could wonder why this idea is still deserving some interest today, when 
the BSSN system is being successfully used in binary-black-holes simulations. Waveform templates 
are currently being extracted for different mass and spin configurations, with an accuracy level 
that depends just on the computational resources (including the use of mesh-refinement and/or 
higher-order finite-difference algorithms). The same is true for neutron stars simulations, where the 
BSSN formalism is currently used for evolving the spacetime geometry [sjj-Q]. But these success 



scenarios have a weak point: the BSSN simulations are basec 



on the combination of the '1+log' 



and 'Gamma-driver' gauge conditions, as proposed in Ref. [j] for the first long-term dynamical 
simulation of a single Black Hole (BH) without excision. 

Concerning BH simulations, we can understand that dealing numerically with collapse singulari- 
ties requires the use of either excision, or time slicing prescriptions with strong singularity-avoidance 
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properties. In the '1+log' case, there is actually a 'limit hypersurface', so that the numerical evolu- 
tion gets safely bounded away from collapse singularities. But singularity- avoidance is a property 
of the time coordinate, which should then be independent of the space coordinates prescription. 
In the spirit of General Relativity, we should expect a gauge-polyvalent numerical code to work as 
well in normal coordinates (zero shift), even if some specific type of time slicing condition (lapse 
choice) is required for BH simulations. Moreover, this requirement should be extended to other 
dynamical choices of the space coordinates. This means that a gauge-polyvalent numerical code 
should also work with alternative shift prescriptions, provided that the proposed choices preserve 
the regularity of the congruence of time lines. And this should be independent of the fact that a 
freezing of the dynamics is obtained or not as a result. These considerations apply 'a fortriori' to 
neutron star simulations without any BH in the final stage, where no singularity is expected to 
form. 

The above proposed gauge-polyvalence requirements, which are in keeping with the spirit of 
General Relativity, may seem too ambitious, allowing for the fact that they are not fulfilled by 
current BH codes. But the need for improvement is even more manifest by looking at the results of 
the gauge- waves test. This test consists in evolving Minkowsky spacetime in non-trivial harmonic 
coordinates, and was devised for cross-comparing the numerical codes performance [8]. In Ref. 
the authors assay different adjustments in order to correct the poor performance of 'plain' BSSN 
codes, which was previously reported in Ref. [9y. They manage to get long-term evolutions for the 
small amplitude case (A = 0.01) with a standard second-order-accurate numerical algorithm. The 
same result was previously achieved by using a fourth-order accurate finite differences scheme [lo| . 
Even in this case, however, the results for the medium amplitude case (A = 0.1) are disappointing. 
More details can be found in a more recent cross-comparison paper [llf], where actually a higher 
benchmark (big amplitude, A = 0.5, devised for testing the non-linear regime) is proposed. 

One could argue that the gauge-waves test is not relevant for real simulations, because periodic 
boundary conditions do not allow constraint violations to propagate out of the computational 
domain 9jj. In BH simulations, however, constraint violations arising inside the horizon can not 
get out, unless all the characteristic speeds of the subsidiary system are adjusted to be greater 
than light speed. As far as this extreme adjustment is not implemented in the current evolution 
formalisms, the gauge- waves test results can be indeed relevant, at least for non-excision BH 
codes. As a result, in keeping with the view expressed in Ref. jlj], we are convinced that either an 
improvement of the current BSSN adjustments or any alternative formulation would be welcome, 
as far as it could contribute to widen the gauge-polyvalence of numerical relativity codes. 
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In this paper we will consider an alternative numerical code consisting in two main ingredients. 
The first one is the Z4 strongly-hyperbolic formulation of the field equations [jjj]. The original 
(second order) version needs no adjustment for the energy and momentum constraints, as far as 
constraint deviations propagate with light speed, although some convenient damping terms have 



been also proposed FJj. We present in Section II a first-order version, which has been adjusted 
for the ordering constraints which arise in the passage from the second-order to the first-order 
formalism. Its flux-Conservative implementation is described in Appendix A. The second ingredient 
is the recently developed FDOC algorithm [la ], which is a (unlimited) finite-difference version 
of the Osher-Chakrabarthy finite- volume algorithm [la ] . along the lines sketched in a previous 
paper [it|]. Although this algorithm, detailed in Appendix B, allows a much higher accuracy, we 
will restrict ourselves here to the simple cases of third and fifth-order accuracy, which have shown 
an outstanding robustness, confirmed by standard tests from Computational Fluid Dynamics, 
including multidimensional shock interactions fis| . 

The results for the gauge-waves test are presented in section where just a small amount 
of dissipation, without any visible dispersion error, shows up after 1000 crossing times, even for 
the high amplitude (A = 0.5) case. Simulations of a 3D BH in normal coordinates are presented 
in section \TV\ where we consider many variants of the 'Bona-Masso' singularity-avoidant prescrip- 
tion 18j j. As expected, the best results for a given resolution are obtained for the choices with a 
limit hypersurface far away from the singularity. For the / = 2/q choice, the BH evolves in normal 
coordinates at least up to 1000 M in a uniform grid with logarithmic space coordinates. This is 
one order of magnitude greater than the normal-coordinates BSSN result, as reported in 7]. 

Concerning the shift conditions, we have tested in Section IVl many explicit first-order prescrip- 
tions in single BH simulations. The idea is just to test the gauge-polyvalence of the code, so no 
physically motivated condition has been imposed, apart from the three-covariance of the shift un- 
der arbitrary time-independent coordinate transformations. Our results confirm that the proposed 
code is not specially tuned for normal coordinates (zero shift). 

II. ADJUSTING THE FIRST-ORDER Z4 FORMALISM 



The Z4 formalism is a covariant extension of the Einstein field Equations, defined as [1 

Rp, + V li Z v + V v Zp = air(T ltt ,--T g^). (1) 



4 



The four vector Z^ is an additional dynamical field, which evolution equations can be obtained 
from (pQ). The solutions of the original Einstein's equations can be recovered when Z^ is a Killing 
vector. In the generic case, the Killing equation has only the trivial solution Z^ = 0, so that true 
Einstein's solutions can be easily recognized. 

The manifestly covariant form ([I]) can be translated into the 3+1 language in the standard way. 
The covariant four-vector Z^ will be decomposed into its space components Zi and the normal 
time component 

e = n lx Z" = a Z° (2) 
where is the unit normal to the t = constant slices. The 3+1 decomposition of (pQ) is given then 



by 
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(d t - C p ) 7y = -2 a Kij (3) 

(d t - Cp) = -Viay + a [Rij + ViZj + 

- 2Kl + (trK-2@) - 8n{Sij -l(trS- r) n 3 } } (4) 

(dt -Cp)@ = | [R + 2 V k Z k + (tr K - 2 6) tr K - tr (K 2 ) - 2 Z k a k /a - Witt] (5) 

(dt - Cp) Zi = a [Vj (Kj - 6i j tr K) + 9,6-2 KJ Z 3 - 9 a,/a - torSi] . (6) 

The evolution system can be completed by providing suitable evolution equations for the lapse 
and shift components. 

d t a = -a 2 Q , dt? = - aQ l (7) 

We will keep open at this point the choice of gauge conditions, so that the gauge-derived quantities 
{Q, Q 1 } can be either a combination of the other dynamical fields or independent quantities with 
their own evolution equation. We are assuming, however, that both lapse and shift are dynamical 
quantities, so that terms involving derivatives of {Q, Q 1 } actually belong to the principal part of 
the evolution system. 

First-order formulation: ordering constraints 

In order to translate the evolution system ([3][7]) into a fully first-order form, the space derivatives 
of the metric components (including lapse and shift) must be introduced as new independent 
quantities: 

Ai = dilna, B k l = d k (i\ D kij = i d k jij . (8) 
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Note that, as far as the new quantities will be computed now through their own evolution equa- 
tions, the original definitions ([8j) must be considered rather as constraints (first-order constraints), 
namely 

A k = A k - d k Ina = (9) 
Bj^Bj-dk? = (10) 
T^kij = D kij - -d k -/ij = 0. (11) 

Note also that we can derive in this way the following set of constraints, related with the ordering 
of second derivatives (ordering constraints): 

Cij = di Aj - dj Ai = di Aj - dj M = , (12) 
Crs* = d r B s i - d s B r l = d r B/ - d s Br* = , (13) 

Crsij = &r *D sij &s T^rij — ®r ^sij B r {j — . (1^0 

The evolution of the lapse and shift space derivatives could be obtained easily, just by taking 
the time derivative of the definitions ([8]) and exchanging the order of time and space derivatives. 
But then the characteristic lines for the transverse-derivative components in (JSj) would be the time 
lines (zero characteristic speed). This can lead to a characteristic degeneracy problem, because 
the characteristic cones of the second-order system dM6|) are basically the light cones JjJ] , and the 
time lines can actually cross the light cones, as it is the case in many black hole simulations. In 
order to avoid this degeneracy problem, we can make use of the shift ordering constraint (I13D for 
obtaining the following evolution equations for the additional quantities ([8]): 

d t A k + di[-(3 l A k + 5 l k (aQ + f3 r A r )\ = B k l A l - trB A k (15) 
dtB k * + dil-tf B k l + S l k (a Q l + p r B r 1 )} = B k l B? - trB B k l (16) 
d t D kij + d^-ftDuj + 8\ {a K {j - 1/2 (By + B^)} ] = B k l D Hj - trB D kij . (17) 

Note that the characteristic lines for the transverse-derivative components are now the normal 
lines (instead of the time lines), so that characteristic crossing is actually avoided. This ordering 
adjustment is crucial for long-term evolution in the dynamical shift case, as it has been yet realized 
in the first-order version of the generalized harmonic formulation . 



Damping terms adjustments 



A further adjustment could be the introduction of some constraint-violation damping terms. 
For the energy-momentum constraints, these terms can be added to the evolution equations ([HE]), 
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as described in Ref. [l4]. 

For the ordering constraints, we can also introduce simple constraint-violation damping terms 
when required. For instance, equation (fT5l) could be modified as follows: 

8 t Ai + d l [-(3 l A l + 5 l t (aQ + (3 r A r )] = B t l A x - trB A i -r ] A l , (18) 

with the damping parameter in the range < r/ <C 1/At. The same pattern could be applied to 
equations (fT6| [T71). 

In order to justify this, let us analyze the resulting evolution equations for the first-order con- 
straints Q. Allowing for (fT5|) . we would get 

d t A k - (? (d r A k - d k Ar) = Bk A r - B/ A k . (19) 

The hyperbolicity of the subsidiary evolution equation (|19p can be analyzed by looking at the 
normal and transverse components of the principal part along any space direction ft, namely 

d t A n -p ± (d n A ± ) = (20) 
d t A ± -(3 n (d n A ± ) = 0, (21) 

with eigenvalues (0, — (3 n ), which is just weakly hyperbolic in the fully degenerate case, that is for 
any space direction orthogonal to the shift vector. Note that this is just the subsidiary system 
governing constraint violations, not the evolution system itself. This means that the main concern 
here is accuracy, rather than stability. But the resulting (linear) secular growth of first-order 
constraint violations may become unacceptable in long-term simulations. 

These considerations explain the importance of adding constraint-damping terms, so that (|15p 
is replaced by (fT8|) . The damping term —r\Ak will appear as a result in the subsidiary system also. 
The linearly growing constraint-violation modes arising from the degenerate coupling in (|20p will 
be kept then under control by these (exponential) damping terms. The same argument applies 
mutatis mutandis to the remaining first-order constraints B^, "D^j . 

Secondary ordering ambiguities 

The shift ordering constraints (|13p can also be used for modifying the first-order version of the 
evolution equation © in the following way 



(d t -Cp) Zi = a[ Vj (Ki j - 5i j trK) + diQ-2 Kf Zj - G At - 8irSi ] - n (dj B£ - $ trB) . (22) 
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Also, the ordering constraints (|14p can be used for selecting a specific first-order form for the three- 
dimensional Ricci tensor appearing in ^ [19[] . This can be any combination of the standard Ricci 
decomposition 

Rij = d^v ij di r kj ~\~ r r k r ^ r r j r ^ (23) 

with the De Donder decomposition 

-Rij = — dk D k ij + 9(j ^j)k ~ ^D r rk Dkij 



which is most commonly used in Numerical Relativity codes. Following Ref. [19J, we will introduce 
an ordering parameter £, so that £ = 1 corresponds to the Ricci decomposition (I23h and £ = — 1 
to the De Donder one (124j) . 

The choices of /i and £ do not affect the characteristic speeds of the evolution system (see 
Appendix A for details), nor the structure of the subsidiary system. In this sense, these are rather 
secondary ordering ambiguities and we will keep these parameters free for the moment, although 
there are some prescriptions that can be theoretically motivated: 

• The choice /i = l/2, £ = —1 allows to recover at the first-order level the equivalence between 
the generalized harmonic formulation and (the second-order version of) the Z4 formalism, 
given by 

Z»= X -Y% a gP° (25) 

(see Appendix A for more details). This can be important, because the harmonic system is 
known to be symmetric hyperbolic. 

• The choice \x = 1 is the only one that ensures the strong hyperbolicity of the Z3 system, 
obtained from the Z4 one by setting 9 = 0. This can be relevant if we are trying to keep 
energy-constraint violations close to zero. Allowing for the quasi-equivalence between the Z3 



and the BSSN systems [191 ]. this adj ustment will affect as well to the first-order version of 
the BSSN system (NOR system [20|) in simulations using dynamical shift conditions. The 



same comment applies to the old 'Bona-Masso' system 



2l|. 



The choice £ = ensures that the first-order version contains only symmetric combinations 
of second derivatives of the space metric. This is a standard symmetrization procedure for 
obtaining a first-order version of a generic second-order equation. 
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In the numerical simulations in this paper, we have taken fj, = 1, £ = — 1, although we have also 
tested other combinations, which also lead to long-term stability. 

III. GAUGE WAVES TEST 

We will begin with a test devised for harmonic coordinates. Let us consider the following line 
element: 

ds 2 = H(x - t){-dt 2 + dx 2 ) + dy 2 + dz 2 , (26) 

where H is an arbitrary function of its argument. One could naively interpret this as the propa- 
gation of an arbitrary wave profile with unit speed. But it is a pure gauge effect, because ([26]) is 
nothing but the Minkowsky metric, written in some non-trivial harmonic coordinates system. 



As proposed in Refs. 
following profile: 



111 ], we will consider the 'gauge waves' line element ([26]) . with the 



H = 1 - A Sin( 2tt(x - t) ) , (27) 



so that the resulting metric is periodic and we can identify for instance the points —0.5 and 0.5 on 
the x axis. This allows to set up periodic boundary conditions in numerical simulations, so that 
the initial profile keeps turning around along the x direction. One can in this way test the long 
term effect of these gauge perturbations. The results show that the linear regime (small amplitude, 
A = 0.01) poses no serious challenge to most Numerical Relativity codes (but see Ref. [l| for the 



BSSN case). Following the recent suggestion in Ref. [Ill ], we will then focus in the medium and 
big amplitude cases (A = 0.1 and A = 0.5, respectively), in order to test the non-linear regime. 
Concerning grid spacing, although Ax = 0.01 would be enough for passing the test in the medium 
amplitude case, the big amplitude one requires more resolution, so we have taken Ax = 0.005 in 
both cases. 

The results of the numerical simulations are displayed in Fig. Q] for the H function (the j xx 
metric component). The left panel shows the medium amplitude case A = 0.1. Only a small 
amount of numerical dissipation is barely visible after 1000 round trips: the third-order-accurate 
finite-difference method gets rid of the dominant dispersion error. For comparison, let us recall 
that the corresponding BSSN simulation crashes before 100 round trips [hJ. The right panel shows 
the same thing for the large amplitude case A = 0.5, well inside the non-linear regime. We see 
some amplitude damping, together with a slight decrease of the mean value of the lapse. 
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FIG. 1: Gauge waves simulation with periodic boundary conditions and sinusoidal initial data for the j xx 
metric component. The resolution is Ax = 0.005 in both cases. The left panel corresponds to the medium 
amplitude case A — 0.1. After 1000 round trips, the evolved profile (cross marks) nearly overlaps the initial 
one (continuous line), which corresponds also with the exact solution. The right panel corresponds to the 
same simulation for the big amplitude case A = 0.5. We see the combination of a slight decrease in the 
mean value plus some amplitude damping. 



Our results are at the same quality level than the ones reported in Ref. ll| for the Flux- 
Conservative generalized-harmonic code Abigail (see also the 'apples with apples' webpage [24]), 
which is remarkable for a test running in strictly harmonic coordinates. We can also com par e 



with the simulations reported in Ref. [23[] for (a specific variant of) the KST evolution system 221 ] . 
Although the gauge wave parametrization is not the standard one, both their 'big amplitude' 
case and their finest resolution are similar to ours. We see a clear phase shift, due to cumulative 
dispersion errors, after about 500 crossing times. We see also a growing amplitude mode, which 
can be moderated with resolution (for the finest one, it just compensates numerical dissipation). 
This can be related with the spurious linear mode that has been reported for harmonic systems 
which are not written in Flux-Conservative form 

We can conclude that there are two specific ingredients in our code that contribute to the 
gauge- wave results in an essential way: the Flux-Conservative form of the equations (see Appendix 
A), which gets rid of the spurious growing amplitude modes, and the third-order accuracy of the 
numerical algorithm, which reduces the dispersion error below the visual detection level in Fig. [H 
even after 1000 crossing times. 
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IV. SINGLE BLACK HOLE TEST: NORMAL COORDINATES 

We will try next to test a Schwarzschild black- hole evolution in normal coordinates (zero shift). 
Harmonic codes are not devised for this gauge choice, so we will compare with BSSN results 
instead. Concerning the time coordinate condition, our choice will be limited by the singularity- 
avoidance requirement, as far as we are not going to excise the black-hole interior. Allowing for 
these considerations, we will determine the gauge evolution equations ([7j) as follows 

Q = f (trK-mQ) , Q* = (0* = 0) , (28) 

where the second gauge parameter m is a feature of the Z4 formalism. We will choose here by 
default m = 2, because the evolution equation for the combination trK— 20, as derived from (jU[5|), 
actually corresponds with the BSSN evolution equation for tr K (see Ref. [19] for the relationship 
between BSSN and Z4 formalisms). 



Concerning the first gauge parameter, we will consider first the '1+log' choice / = 2/a [25f] , 
which is the one used in current binary BH simulations in the BSSN formalism. The name comes 
from the resulting form of the lapse, after integrating the evolution equation (O [7]) with the 
prescription (|28p for true Einstein's solutions (O = 0): 

a = a + In (7/70) , (29) 

where is the space volume element. It follows from ([29]) that the coordinate time evolution 
stops at some limit hypersurface, before even getting close to the collapse singularity. This happens 
when 

VtTto = exp(-a /2) , (30) 

that is well before the vanishing of the space volume element: the initial lapse value is usually close 
to one, so that the final volume element is still about a 60% of the initial one. This can explain 
the robustness of the 1+log choice in current black-hole simulations. 

We will consider as usual initial data on a time-symmetric time slice (Kij = 0) with the intrinsic 
metric given in isotropic coordinates: 

fa = (! + 1:) 4 % • ( 31 ) 

This is the usual 'puncture' metric, with the apparent horizon at r = m/2: the interior region is 
isometric to the exterior one, so that the r = singularity is actually the image of space infinity. 
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We prefer, however, to deal with non-singular initial data. We will then replace the constant mass 
profile in interior region r < M/2 by a suitable profile m(r), so that the interior metric corresponds 
to a scalar field matter content. Of course, the scalar field itself must be evolved consistently there 
(see Appendix C for details). A previous implementation of the same idea, with dust interior 
metrics, can be found in Ref. [261 ] . 
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FIG. 2: Plots of the lapse profiles at t = 20M and t — 40 M. The results for the third-order accurate 
algorithm (continuous lines) are compared with those for the fifth-order algorithm (dotted lines) for the 
same resolution (h — 0.1M). We have also included for comparison one extra line, corresponding to the 
third-order results with h — 0.05 M, computed in a reduced mesh. Increasing resolution leads to a slope 
steepening and a slower propagation of the collapse front. In this sense, as we can see for t = 20 M, switching 
to the fifth-order algorithm while keeping h = 0.1 M amounts to doubling the resolution for the third-order 
algorithm. 

We have performed a numerical simulation for the / = 2/ a case with a uniform grid with 
resolution h = 0.1 M, extending up to r = 20 M (no mesh-refinement). We have used the third and 
fifth-order FDOC algorithms, as described in Appendix B, with the optimal dissipation parameters 
for each case. The results for the lapse profile are shown in Fig. [2] at t = 20 M an t = 40 M. We see 
in both cases that the higher order algorithm leads to steeper profiles and a slower propagation of 
the collapse front. Note that the differences in the front propagation speed keep growing in time, 
although the third-order plot at t = 40 M is clearly affected by the vicinity of the outer boundary. 
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This fact does not affect the code stability, as far as we can proceed with the simulations beyond 
t = 50 M, when the collapse front gets out of the computational domain (beyond t = 60 M in 
the higher-order simulations). Note that the corresponding BSSN simulations (/ = 2/a in normal 
coordinates) are reported to crash at about t = 40 M [7]. 

We have added for comparison an extra plot in Fig. [2j with the results at t = 20 M of a third- 
order simulation with double resolution {h = 0.05 M), obtained in a smaller computational domain 
(extending up to 10 M). Both the position and the slope of the collapse front coincide with those 
of the fitfth-order algorithm with h = 0.1 M. In this case, switching to the higher-order algorithm 
amounts to doubling accuracy. Note, however, that higher-order algorithms are known to be less 
robust [3]. Moreover, as the profiles steepen, the risk of under-resolution at the collapse front 
increases. We have found that a fifth-order algorithm is a convenient trade-off for our h = 0.1M 
resolution in isotropic coordinates. 

We have also explored other slicing prescriptions with limit surfaces closer to the singularity, 
as described in Table HI Note that in these cases the collapse front gets steeper than the one 
shown in Fig. [2] for the standard / = 2/a case with the same resolution. This poses an extra 
challenge to numerical algorithms, so we have switched to the third-order-accurate one for the sake 
of robustness. In all cases, the simulations reached t = 50 M without problem, meaning that the 
collapse front has get out of the computational domain. It follows that the standard prescription 
f = 2/a, although it leads actually to smoother profiles, is not crucial for code stability. 



f 


2/a 


1+1/a 


1/2+1/a 


l/a 


V7/70 


61% 


50% 


44% 


37% 



TABLE I: Different prescriptions for the gauge parameter /, with the corresponding values of the residual 
volume element at the limit surface (normal coordinates), assuming a unit value of the initial lapse. 



The results shown in Fig. [2] compare with the ones in Ref. [27J], obtained with (a second-order 
version of) the old Bona-Masso formalism. We see the same kind of steep profiles, produced by the 



well known slice-stretching mechanism 28]. This poses a challenge to standard numerical methods: 
in Ref. 27] F inite- Volume methods where used, including slope limiters. Our FDOC algorithm 
(see Ref. 15J for details) can also be interpreted as an efficient Finite-Differences (unlimited) 
version of the Osher-Chakrabarthy Finite- Volume algorithm [16]. Note however that in Ref. [27]. 
like in the BSSN case, a conformal decomposition of the space metric was considered, and an 
spurious (numerical) trace mode arise in the trace-free part of the extrinsic curvature. An additional 
mechanism for resetting this trace to zero was actually required for stability. In our (first-order) 
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Z4 simulations, both the plain space metric and extrinsic curvature can be used directly instead, 
without requiring any such trace-cleaning mechanisms. 

Let us take one further step. Note that the lifetime of our isotropic coordinates simulations 
(with no shift) is clearly limited by the vicinity of the boundary (at r = 20 M). At this point, we 
can appeal to space coordinates freedom, switching to some logarithmic coordinates, as defined by 

R = L sinh{r/L) , (32) 

where R is the new radial coordinate and L some length scale factor. This configuration suggests 
using the third-order algorithm because of its higher robustness. We have performed a long-term 
numerical simulation for the / = 2/ a case, with L = 1.5 M, so that R = 2QM in these logarithmic 
coordinates corresponds to about r = 463.000 M in the original isotropic coordinates. In this way, 
as shown in Fig. [3l the collapse front is safely away from the boundary, even at very late times. We 
stopped our code at t = 1000 M, without any sign of instability. This provides a new benchmark 
for Numerical Relativity codes: a long-term simulation of a single black-hole, without excision, in 
normal coordinates (zero shift). Moreover, it shows that a non-trivial shift prescription is not a 
requisite for code stability in BH simulations. 



V. SINGLE BLACK HOLE TEST: FIRST-ORDER SHIFT CONDITIONS 

Looking at the results of the previous Section, one can wonder wether our code is just tuned for 
normal coordinates. This is why we will consider here again BH simulations, but this time with 
some non-trivial shift prescriptions. The idea is just to test some simple cases in order to show 
the gauge-poly valence of the code. For the sake of simplicity, we will consider here just first order 
shift prescriptions, meaning that the source terms (Q, Q l ) in the gauge evolutions (|T|) are algebraic 
combinations of the remaining dynamical fields. To be more specific, we shall keep considering 
slicing conditions defined by 

Q = -f3 k /a A k + f (trK - m&) , (33) 



together with dynamical shift prescriptions, defined by different choices of Q l 



We will 



First-order shift prescriptions have been yet considered at the theoretical level 
introduce here an additional requirement, which follows when realizing that, allowing for the 3+1 
decomposition of the line element 

ds 2 = -a 2 dt 2 + 7y (dx i + ftdt) (dx j + fidt) , (34) 
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FIG. 3: Plot of the lapse function for a single BH at t = 1000 M in normal coordinates. Only one of every ten 
points is shown along each direction. The third-order accurate algorithm has been used with (3 = 1/12 and 
a space resolution h — 0.1 M. The profile is steep, but smooth: no sign of instability appears. Small riddles, 
barely visible on the top of the collapse front, signal some lack of resolution because of the logarithmic 
character of the grid. The dynamical zone is safely away from the boundaries. 

the shift behaves as a vector under (time independent) transformations of the space coordinates. 
We will impose then that its evolution equation, and then Q % , is also three-covariant. 

This three-covariance requirement could seem a trivial one. But note that the harmonic shift 
conditions, derived from 



are not three-covariant (the box here stands for the wave operator acting on scalars). In the 3+1 
language, (f35l) can be translated as 



where the non-covariance comes from the space-derivatives terms. 

Concerning the advection term, a three-covariant alternative would be provided either by the 
Lie-derivative term 



□ x l = 



(35) 



(36) 



(37) 
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or by the three-covariant derivative term 

P k V k (f3 l /a) = l/a[p k B k * - (3 l (3 k A k + T) k (3^ k ] . (38) 

We have tested both cases in our numerical simulations. 

Concerning the last term in (|36p , we can take any combination of A 1 , Z l and the vectors obtained 
form the space metric derivatives after subtracting their initial values, namely: 

Di — Di | t= o , Ei — Ei |t=o • (39) 

This is because the additional terms arising in the transformation of the non-covariant quantities 
(D{, Ei) depend only on the space coordinates transformation, which is assumed to be time- 
independent. Note that, for the conformal contracted-Gamma combination 

Ti = 2Ei-^Di , (40) 



the subtracted terms actually vanish in simulations starting from the isotropic initial metric ([31 
Of course, the same remark applies to the BSSN Gamma quantity, namely [19 ] 



T i + 2Z l . (41) 



We have considered the following combinations: 



• a 2 

SI: d t (3 l = — A 1 - a Q (i l (42) 



o 



2 



S2: d t if = ^- A t + P k B h i + T A jk pp k -aQp t (43) 



2 



n 



2 



53: d t p l = — r + p k B k * + T) k p{3 k -aQP\ (44) 

where SI corresponds to the Lie-derivative term (|37p and the remaining two choices to the covariant 
advection term (|38p . with different combinations of the first-order vector fields. 

We have obtained stable evolution in all cases, with the simulations lasting up to the point 
when the collapse front crosses the outer boundary (about t = 50 M). We can see in Fig. [4] the 
lapse and shift profiles in the SI and the S3 cases (S2 is very similar to SI). The shift profiles 
are modulated by the lapse ones, so that the shift goes to zero in the collapsed regions. This is a 
consequence of the term —a Q f3 % in the shift evolution equation, devised for getting finite values 
of the combination (3 l /a. In the non-collapsed region, SI leads to a higher shift profile, which 
spreads out with time, whereas S3 leads to a lower profile, which starts diminishing after the initial 
growing. Allowing for (144j) . this indicates that the conformal gamma quantity Tj is driven to zero. 
The lapse slopes are also slightly softened in the S3 case. 
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FIG. 4: Plot of the lapse and shift profiles at t = 20 M (continuous lines) and t = 40 M (dotted lines). The 
plots are shown along the main diagonal of the computational domain, in order to keep the outer boundary 
out of the dynamical zone. In the SI case (left panel), after the initial growing, the maximum shift value 
keeps constant. In the S3 case (right panel), it clearly diminishes with time. 

These results confirm that the code stability is not linked to any particular shift prescription, 
as we can combine different source terms in the shift evolution equation, leading to different lapse 
and shift profiles. 

VI. CONCLUSIONS AND OUTLOOK 



We have shown in this paper how a first-order flux-conservative version of the Z4 formalism can 
be adjusted for dealing with the ordering constraints, and then implemented in a numerical code 
by means of a robust, cost-efficient, finite-difference formula. The resulting scheme has been tested 
in a demanding harmonic-coordinates scenario: the gauge-waves testbed. The code performance 



111 ], even in the highly non- linear 



compares well with the best harmonic-code results for this test 
regime (50% amplitude case). This is in contrast with the well-known problems of BSSN-based 

n n 

codes with the gauge- waves test [I] [8j]. 

The code has also been tested in non-excision BH evolutions, where singularity-avoidance is a 
requirement. Our results confirm the robustness of the code for many different choices of dynamical 
lapse and shift prescriptions. In the normal coordinates case (zero shift), our results set up a new 
benchmark, by evolving the BH up to 1000 M without any sign of instability. This improves 
the reported BSSN result by one order of magnitude (Harmonic codes are not devised for normal 
coordinates). More important, this shows that a specific shift choice is not crucial for code stability, 
even in non-excision BH simulations. This is confirmed by our shift simulations, where different 



17 



covariant evolution equations for the shift lead also to stable numerical evolution. 

In spite of the encouraging performance in these basic tests, we still are on the way towards 
a gauge-polyvalent code, as pointed out by the title of this paper. More technical developments 
on the numerical part are required: mesh refinement, improved boundary treatment, etc. On the 
theoretical side, as far as the shift prescription is no longer determined by numerical stability, we 
can explore shift choices from the physical point of view, adapting our space coordinates system 
to the features of every particular problem. We are currently working in these directions. 
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Appendix A: Flux-Conservative evolution equations 

We will write the first-order evolution system in a balance-law form. For a generic quantity u, 
this leads to 

dt u + d k F\u) = S{u) , (A.l) 

where the Flux F k (u) and Source terms S(u) can depend on the full set of dynamical fields in an 
algebraic way. In the case of the space-derivatives fields, their evolution equations (|15lll7p are yet 
in the balance-law form (jA.ip . Note however that any damping terms of the form described in (|18|) 
will contribute both to the Flux and the Source terms in a simple way. 
The metric evolution equation (|3|) will be written in the form 

dt la = 2 (3 k D kij + Bij + B ji -2a Ky , (A.2) 

so that it is free of any Flux terms. The remaining (non-trivial) evolution equations (JU- [6|) require 



18 



a more detailed development. We will expand first the Flux terms in the following way: 

dtKij + d k [-(3 k Kij + a X k l3 } = S(K tJ ) (A.3) 

dtZi + d k [-(3 k Zi + a {-K k i + S k i(trK - 9)} (A.4) 

+fi (Bi k - 5 t k trB) } = S{Zi) 

d t & + d k [-(3 k e + a(D k - E k - Z k )]= S(Q) (A.5) 

where we have used the shortcuts Di = D ik k and Ei = D k ki , and 

\\, = Eft* - \ (1 + (A/ + D 3 k ) + \ 5 k t [A, + Dj - (1 - Ej -2Zj] (A.6) 

+ X -b k 3 [At + Di-(l-g) Ei-2 Zi]. 

The Source terms S(u) do not belong to the principal part and will be displayed later. Let us 
focus for the moment in the hyperbolicity analysis, by selecting a specific space direction ft, so that 
the corresponding characteristic matrix is 

d F n 

A n = ^—, (A.7) 

where the symbol n replacing an index stands for the projection along the selected direction n. We 
can get by inspection the following (partial) set of eigenfields, independently of the gauge choice: 

• Transverse derivatives: 

A ± , B ± \ D±ij , (A.8) 

propagating along the normal lines (characteristic speed —f3 n ). The symbol _L replacing an 
index means the projection orthogonal to n. 

• Light-cone eigenfields, given by the pairs 

F n [D n±± ] ± F n [K ±± ] (A.9) 
-F n [Z ± ] ± F n [K n± ] (A.10) 
F n [D n -E n -Z n ] ± F n [Q] (A.ll) 

with characteristic speed —j3 n ± a , respectively. 

Note that the eigenvector expressions given above, in terms of the Fluxes, are valid for any choice 
of the ordering parameters fi and £. Only the detailed expression of the eigenvectors, obtained from 
the Flux definitions, is affected by these parameter choices. For instance 

F n [D n — E n — Z n ] = —(3 n [D n — E n — Z n ] + a + (n — l)tr(B ±± ) . (A.12) 
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Any value /i / 1 implies that the characteristic matrix of the Z3 system, obtained by removing the 
variable 9 from our Z4 evolution system [191, can not be fully diagonalized in the dynamical shift 
case. Of course, the hyperbolicity analysis can not be completed until we get suitable coordinate 
conditions, amounting to some prescription for the lapse and shift sources Q and Qi, respectively. 
But the subset of eigenvectors given here is gauge independent: non-diagonal blocs can not be 
fixed a posteriori by the coordinates choice. 

The detailed expressions for the eigenvectors can be relevant when trying to compare with 
related formulations. For instance, a straightforward calculation shows that the eigenvectors (|A.9l 
lA.llj) can be matched to the corresponding ones in the harmonic formalism if and only if 

£ = -1, /i = l/2. (A.13) 

This shows that different requirements can point to different choices of these ordering parameters. 
We prefer then to leave this choice open for future applications. Concerning the simulations in this 
paper, we have taken £ = — 1 , fj, = 1 . 

Finally, we give for completeness the Source terms, namely: 

S(K l3 ) = -K id trB + K ik B k + K jk B* + a (1 + £) [-A k Y k i3 + h^A* D 3 + A 3 A)] 
+ i (1 - £) [A k D k i3 - X -{A 3 (2 Ei - A) + A (2 E 3 - D 3 )} 
+ 2 (Di r m D r m j + Dj r m D 7 ' m i) — 2 E k (Di 3 k + D 3 i k )] 
+ (D k + A k - 2 Z k ) T k i 3 - T k m3 T m ki - (Ai Z 3 + A 3 Z£) - 2 K k i K kj 
+ (trK - 2 G) K i3 } - 8 vr a [S l3 - ± (trS - r) li3 ] (A.14) 
S(Zi) = -Zi trB + Z k B { k + a [A, (trK - 2 9) - A k K k i - K k r T r ki + K k \ (D k - 2 Z k )\ 

-8 vr a ^ (A.15) 
S(0) = -9 trB + |[2 A k (D k - E k - 2 Z k ) + D k rs T k rs - D k (D k - 2 Z k ) - K k r K r k 

+ trK (trK - 2 9)] - 8 vr a r . (A.16) 



Appendix B: Finite-differences implementation 



We follow the well-known method-of-lines (MoL [30j] ) in order to deal separately with the space 
and the time discretization. Concerning the time discretization, we use the following third-order- 
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accurate Runge-Kutta algorithm 



u* = u n + At rhs( u n ) 

u** = ^u n + - [u* + At rhs( u* )] (B.l) 
u n+l = \u n + ~ [u** + At rhs{ u** )] , 



which is strong-stability-preserving (SSP 31]), where we have used as a shorthand 

rhs{ u ) = -d k F k (u) + S(u) . (B.2) 

The Flux derivatives appearing in (|B.2|) will be discretized by using the finite-difference formula 
proposed in Ref. [15] (FDOC algorithm). For instance the derivative of F x {u) will be represented 
as 

3 x Ff = C 2m Ff + (-l) m /?(Ax) 2m D™D m " 1 (A J _ 1/2 D^ Uj ) , (B.3) 

where C 2m is the 2mth-order-accurate central difference operator and D± are the standard finite 
difference operators. We have also noted 

Xj-i/2 = max(X j ,X j ^ l ) , (B.4) 

where Xj stands here for the local characteristic radius (the highest characteristic speed, typically 
the gauge speed). 

Note that the second term in the finite-difference formula (|B.3P is actually a dissipation operator 
of order 2m acting on (A u), so it could be regarded at the first sight as a mere generalization of the 
standard Kreiss-Oliger artificial viscosity operators [32]. This is not the case: the formula (|B.3j) 
can be instead derived in a finite-volume framework, when combining the local-Lax-Friedrichs flux 
H with the (—ed) Osher-Chataha^ 9uX integration Q (see Ref. Q for 



formula 



details, including the optimal values of the (5 parameter). 

Note that, contrary to the standard Kreiss-Oliger approach, the dissipation term is such that 
the accuracy of the first (centered derivatives) term in (|B.3P is reduced by one order: the resulting 
FDOC algorithm accuracy is always of an odd order. This is important for code robustness. The 
algorithms (IB.3D can be shown to keep monotonicity even for remarkably high compression factors 
(defined as the ratio between two neighbor slopes along a given direction) which is what is 
actually required in view of the steep profiles shown for instance in Fig. [2j 

The space accuracy of the scheme (IB. 3ft is 2m — 1, with an stencil of 2m+l points. We have used 
in this paper both the third-order and fifth-order accurate methods, for which the optimal values 
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of the dissipation parameter are (3 = 1/12, j3 = 2/75, respectively [151 ] . In the fifth order case, 
we have a seven-point stencil and the dissipation term corresponds to a sixth derivative, as in the 
advanced finite-difference schemes used in Ref. 



341 ], The robustness of the proposed algorithms, 



with compression factors of 5 and 3 respectively, makes them very convenient for steep-gradient 
scenarios, such us the ones arising in black-hole simulations, where slice-stretching threatens the 
stability of more standard finite-difference algorithms [281 ] . 

No sophisticated numerical tools (mesh refinement, algorithm-switching for the advection terms, 
etc) have been incorporated to our code at this point, when we are facing just test simulations. 
Concerning the boundary treatment, we simply choose at the points next to the boundary the most 
accurate centered algorithm compatible with the available stencil there. When it comes to the last 
point, we can either copy the neighbor value or propagate it out with the maximum propagation 
speed (by means of a ID advection equation). The idea is to keep the numerical code as simple as 
possible in order to test here just the basic algorithm in a clean way. 



Appendix C: Scalar field stuffing 

Let us consider the stress-energy tensor 

T ab = $a$b- 1/2 (g Cd $c $d) 9ab , (C.l) 

where we have noted $ a = <9 a <3?, corresponding to a scalar field matter content. The 3+1 decom- 
position of (jC.ip is given by 

r= l/2($ n 2 + 7 fc ^ fc ^), # = % = + 1/2 (*„ 2 - t* 1 ** 7»; , (C.2) 

where $ n stands for the normal time derivative: 

(d t - (5 k d k ) $ = (C.3) 

The quantities (|C.2p appear as source terms in the field equations dMBJ). 

The stress-energy conservation amounts to the evolution equation for the scalar field, which is 
just the scalar wave equation. In the 3+1 language, it translates into the Flux-conservative form: 

dt Wl$n] + d h [^(-/3 fc a> n + a 7 fc %)] =0. (C.4) 

A fully first-order system may be obtained by considering the space derivatives ^ as independent 
dynamical fields, as we did for the metric space derivatives. 
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Concerning the initial data, we must solve the energy- momentum constraints. They can be 
obtained by setting both and Zi to zero in (|5j E|). In the time-symmetric case (Kij = 0), this 
amounts to 

R=16ttt, Si = <f> n <S>i = Q. (C.5) 

The momentum constraint will be satisfied by taking $ (and then $j) to be zero everywhere on 
the initial time slice. Concerning the energy constraint, we will consider the line element (|3ip with 
m = m(r). We assume a constant mass value m = M for the black-hole exterior, so that the 
energy constraint in (|C,5P will be satisfied with r = there. 

In the interior region, the energy constraint will translate instead into the equation 

m" = -2vrr (<D n ) 2 (1 + ^) 5 , (C.6) 

which can be interpreted as providing the initial $ n value for any convex (to" < 0) mass profile. 
Of course, some regularity conditions both at the center and at the matching point r$ must be 
assumed. Allowing for (|C.6p . we have taken 

to = to" = (r = 0) 
to = M, m! = to" = (r = tq) . 

Note that, allowing for (|C.6h . these matching conditions ensure just the continuity of $ n , not 
its smoothness. This can cause some numerical error, as we are currently evolving $ n through 
the differential equation (|C.4p . If this is a problem, we can demand the vanishing of additional 
derivatives of the mass function m(r), both at the origin and at the matching point (this is actually 
the case in our shift simulations). This is not required in the standard case (/ = 2/a, normal 
coordinates), where we have used a simple profile, with the matching point at the apparent horizon 
(r = M/2), given by 

m(r) = 4r - 4/M[r 2 + (M/2vr) 2 sm 2 (2vrr/M) ] . (C.7) 
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